% BL_OutpuGap_REStat_1 - Figures 1, 4, 5, and 7 and Table 2
%
% Replication files for the paper "Measuring the output gap with large datasets"
% by Matteo Barigozzi and Matteo Luciani, forthcoming on the Review of economics and statistics.
% 

% Written by Matteo Luciani, matteoluciani@yahoo.it

clear all; close all; clc; 

ML_graph_options


STIMA=0;  % if STIMA==1 it estimates the whole model including confidence bands
SALVA=0;  % if SALVA==1 it saves the results
STAMPA=0;  % if STAMPA==1 it prints the charts

if STIMA==0 
   

elseif STIMA==1  % -------------------------------------------------------- % Estimate everything
 
    %%% ============================================= %%%
    %%%  Set-up settings for benchmark specification  %%%
    %%% ============================================= %%%
    Setting.giorno='20190329';                                              % RESTAT SUBMISSION
    Setting.gm=datenum(str2num(Setting.giorno(1:4)),...                     % -----------------
        str2num(Setting.giorno(5:6)),str2num(Setting.giorno(7:8)));         % -----------------
    Setting.Filename=['USDB_Haver_' num2str(Setting.gm)];                   % file name    
    Setting.tresh=10^(-2);                                                  % tolerance treshold for EM algorithm
    Setting.star=10^(-5);                                                   % Initial variance R
    Setting.maxiter=50;                                                     % max number iteratioin EM algorithm
    Setting.trans=3;                                                        % data transformation   
    Setting.out=2;                                                          % outlier treatment
    Setting.q=4;                                                            % number of factors
    Setting.s=1;                                                            % lags in the factor loadings
    Setting.d=Setting.q-1;                                                  % number of common cycles
    Setting.p=3;                                                            % lags VAR
    Setting.det=1;                                                          % parameters DFM
    Setting.GDO=1;                                                          % impose GDO restrictions
    Setting.model='VAR';                                                    % determines law of motion for the factors
    Setting.m=0;                                                            % parameters Robinson-Yao-Zhang
    Setting.cc=10;                                                          % obs to exclude because of initial condition
    Setting.TR1=[6 15 72 73 76:80 87];                                      % Variables for which I overwrite the trend test
    Setting.I0=[1:2 36 38 40 43 68:71 75];                                  % restrictions for EM algorithm
    Setting.I1=[];                                                          % -----------------------------
    Setting.rr=ones(Setting.q,1);                                           % -----------------------------
    Setting.TV.id={1:2,75};                                                 % time varying parameters
    Setting.TV.Type={['trend';'none '],'mean'};                             % -----------------------
    Setting.TV.q0=[10^(-3), 10^(-2)];                                       % initial variance for TV states
    Setting.nboot=1000;                                                     % Number of bootstrap
    
    save RestatSetting Setting
    
    %%% =========== %%%
    %%%  Load Data  %%%
    %%% =========== %%%
    [Y, Label, Name, Dates, cd, CBO, idcat, NameC] = ...                    % Load data in levels
        ML_ReadDataHaver(Filename,trans,out,cut,StartYear);                 % -------------------
    
    %%% =============================== %%%
    %%%  Initialize the model with PCA  %%%
    %%% =============================== %%%
    y=ML_diff(Y);                                                           % Data in 1st Differences
    [T,N]=size(y);                                                          % size of the panel
    [yy, my, sy]=ML_Standardize(y);                                         % Standardize
    [f,lambda]=ML_efactors2(yy,q,2);                                        % estimate factor loadings, aka DFM on \Delta y_t
    TT=(1:T+1)';                                                            % time trend
    X=NaN(T+1,N); bt=X; b=zeros(N,2);                                       % preallocates variables for detrending
    J=BL_TestLinearTrend(y); J(TR1)=1;                                      % Identify variables to be detrended
    [X(:,J),bt(:,J),b(J,:)]=ML_detrend(Y(:,J));                             % Detrend variables to be detrended
    X(:,~J)=Y(:,~J)-repmat(mean(Y(:,~J)),T+1,1);                            % Demean variables not to be detrended
    bt(:,~J)=repmat(mean(Y(:,~J)),T+1,1);                                   % ---------------------
    b(~J,1)=mean(Y(:,~J))';b(~J,2)=0;                                       % ---------------------
    if GDO==1 % ----------------------------------------------------------- % Restrictions for GDO
        lambda(1:2,:)=repmat(mean(lambda(1:2,:)),2,1);                      % same loadings
        b(1:2,2)=mean(b(1:2,2));                                            % same slope, different constant
        bt(:,1:2)=repmat(b(1:2,1)',T+1,1)+TT*b(1:2,2)';                     % same linear trend
        X(:,1:2)=(Y(:,1:2)-bt(:,1:2));                                      % Detrended GDP GDI
        sy(1:2)=repmat(mean(sy(1:2)),1,2);                                  % same std for standardization
    end        % -------------------------------------------------------------- %
    Z=X./repmat(sy,size(X,1),1);                                            % BLL - Standardization
    F=Z*lambda/N;                                                           % Factors in levels as in BLL
    [A0,v,AL]=ML_VAR(F,p,1);                                                % Estimate VAR
    xi=Z-F*lambda';                                                         % idiosyncratic component
    
    %%% ============================ %%%
    %%%  Estimate the model with EM  %%%
    %%% ============================ %%%
    Z2=(Y-repmat(b(:,1)',size(Y,1),1))./repmat(sy,size(Y,1),1);             % BL standardizations: (Y_t - a), i.e. approx. centered levels, divided by sy
    b1=b./repmat(sy',1,2); b1(:,1)=0;                                       % divide slope by sy, constant=0;
    NSDFM_SS=ML_NSDFM_SS_GDO_TV(AL,lambda,v,q,s,xi,F,rr,I0,I1,b1,A0,TV,star);   % State-space representation
    EM=ML_NSDFM_EM_GDO_TV(Z2,NSDFM_SS,iter,tresh,cc,GDO);                   % EM-Algorithm
    T2=size(Z2,1); start=1;
    BT=TT*EM.b2';                                                           % ML estimates of linera trend
    isTV=find(NSDFM_SS.type2==10)+max(p*q,q*(s+1));                         % identifies TV coefficients
    for ii=1:length(isTV) % --------------------------------------------------- % TV slopes or means
        id=find(EM.Lambda(:,isTV(ii)));
        BT(:,id)=EM.xitT(:,isTV(ii))*ones(1,length(id));
    end                   % --------------------------------------------------- %
    FF=EM.xitT(:,1:q); FF1=ML_lag2(FF,s,0);                                 % Store ML Factors
    L=EM.Lambda(:,1:(s+1)*q);                                               % Store ML loadings
    L2=NaN(N,q,s+1); for ss=1:s+1; L2(:,:,ss)=L(:,(ss-1)*q+1:ss*q); end     % -----------------
    
    
    %%% =================== %%%
    %%%  Common Components  %%%
    %%% =================== %%%
    T3=length(FF1(cc:end,:));
    SY=repmat(sy,T3,1);
    MY=repmat(b(:,1)',T3,1);
    start2=start+s+cc-1;
    Y2=Y(start2:end,:);
    chi=(FF1(cc:end,:)*L'+BT(start2:end,:)).*SY+MY;                         % common component
    zeta=Y2-chi;                                                            % idiosyncratic component
    
    %%% ================================= %%%
    %%%  Estimate trend with PCA on LRCV  %%%
    %%% ================================= %%%
    [V, D]= eig(cov(FF(cc:end,:)));                                         % eigenvalue and eigenvenctors of VCV
    [~, t2]=sort(diag(D),'descend'); V=V(:,t2); D=D(t2,t2);                 % sort eigenvalues and eigenvectors
    FFt = FF(cc:end,:)*V(:,1:q-d);                                          % common trend
    FFc = FF(cc:end,:)*V(:,q-d+1:q);                                        % common cycle
    FFt1=[]; FFc1=[];
    for ss=1:s+1
        FFt1=cat(2,FFt1,FFt(s-ss+2:end-ss+1,:));                            % lagged common trend
        FFc1=cat(2,FFc1,FFc(s-ss+2:end-ss+1,:));                            % lagged common cycle
        lambdat(:,(ss-1)*(q-d)+1:ss*(q-d)) = L2(:,:,ss)*V(:,1:q-d);         % trend loadings
        lambdac(:,(ss-1)*d+1:ss*d) = L2(:,:,ss)*V(:,q-d+1:q);               % cycle loadings
        LcT(:,:,ss)=L2(:,:,ss)*V(:,q-d+1:q);
    end
    
    chit = (FFt1*lambdat'+BT(start2:end,:)).*SY+MY;                         % common stochastic trend plus linear trend
    chist = (FFt1*lambdat').*SY+MY;                                         % common stochastic trend
    chic = (FFc1*lambdac').*SY;                                             % common temporary
    chint = (FFt1*lambdat').*SY;                                            % common stochastic trend
    chilt = (BT(start2:end,:)).*SY;                                         % linear trend
    
    %%% ============================ %%%
    %%%  Bootstrap confidence bands  %%%
    %%% ============================ %%%
    eta=Z2-EM.xitT*EM.Lambda';
    chiB=NaN*repmat(chi,1,1,nboot); chicB=chiB; chitB=chiB;
    tic
    type2=NSDFM_SS.type2;
    
    parfor bb=1:nboot; disp(bb)
        [Z2s,xitTs]=ML_NSDFM_DGP_GDO_TV(T+1,N,q,s,p,EM,eta,type2,[1 0],Block);
        Boot=ML_NSDFM_Boot_GDO_TV(Z2s,NSDFM_SS,GDO,iter,tresh,cc,EM,xitTs,d);
        chiB(:,:,bb)=Boot.chi.*SY+MY;
        chicB(:,:,bb)=Boot.chic.*SY;
        chitB(:,:,bb)=Boot.chit.*SY+MY;
    end
    toc
    
    %%% ============= %%%
    %%% CBO Estimates %%%
    %%% ============= %%%
    J=ismember(CBO(:,1),Dates);
    cbot=100*log(CBO(J,2));                                                     % Potential Output CBO
    cbost=cbot-bt(:,1);                                                         % Potential Output CBO
    cboc=Y(:,1)-cbot;                                                           % Output gap CBO
    cbonr=CBO(J,3);                                                             % Natural Rate CBO
    cboug=Y(:,75)-cbonr;                                                        % Unemployment gap CBO        
    
    %%% ====================== %%%
    %%% == ---------------- == %%%
    %%% == Plotting Results == %%%
    %%% == ---------------- == %%%
    %%% ====================== %%%
    
    alpha1=16; alpha2=32;
    Dates2=Dates(start2:end); Dates2d=Dates2(2:end); Dates4q=Dates2(5:end);     % New Dates for plots
    j0=find(year(Dates2)==1970,1);                                              % Strarting point for all graphs
    colore={[0 0.45 0.74],[0.64 0.08 0.18],[0.85 0.33 0.1],...                  % colors for gaps blue, red, orange
        [0.93 0.69 0.13],[0 0 0],[0.47 0.67 0.19],[0.49 0.18 0.56]};            % yellow,black, green, purple
    nm=size(chic,2);                                                            % number of possible decompositions
    qsp=[num2str(q) num2str(s) num2str(p)];
    DD=Dates2(j0:end);

    %%% ====================== %%%
    %%% == ---------------- == %%%
    %%% ==   SAVE RESULTS   == %%%
    %%% == ---------------- == %%%
    %%% ====================== %%%
    
    if SALVA==1
        tic
        save RestatBenchmarkALL.mat         % large file containing the whole distribution
        
        RESTAT.chi=chi;         RESTAT.chic=chic;   RESTAT.chit=chit;
        RESTAT.chint = chint;   RESTAT.chilt = chilt;   RESTAT.zeta=zeta;
        RESTAT.ogapB=squeeze(chicB(:,1,:));     RESTAT.potB=squeeze(chitB(:,1,:));
        RESTAT.ugapB=squeeze(chicB(:,75,:));    RESTAT.nrB=squeeze(chitB(:,75,:));
        RESTAT.alpha=[alpha1 alpha2];
        RESTAT.FF=FF;           RESTAT.L=L;         RESTAT.BT=BT;
        RESTAT.FFt=FFt;         RESTAT.lambdat=lambdat;
        RESTAT.FFc=FFc;         RESTAT.lambdac=lambdac;               
        save RestatBenchmark2 RESTAT;   % smaller file containing just selected percentiles
        toc
    end

end


%%% ========= %%%
%%%  Table 2  %%%
%%% ========= %%%
ID=[1 40 43 41 69 75 72 73 82 97 103];
Tabella=[lambdat(ID,:) lambdac(ID,[3 6 1 4]) -lambdac(ID,[2 5])];

LS={'k-','k--','k:','k-.'};
nomefile='Restat_00bw_';


lg1={'CBO','BL'};
lg2={'GDP','CBO','BL'};



%%% ======================== %%%
%%%  Figure 1 & 5 left plot  %%%
%%% ======================== %%%
ZZ=[cboc(start2:end) chic(:,1)]; ZZ=ZZ(j0:end,:);                           % Output Gap - levels
ZZb{1}=BL_Band(chic(j0:end,1),squeeze(chicB(j0:end,1,:)),alpha1);           % -------------------
ZZb{2}=BL_Band(chic(j0:end,1),squeeze(chicB(j0:end,1,:)),alpha2);           % -------------------
SizeXY=ML_SizeXY(DD,ML_minmax([ZZ ZZb{1}]),1,5,[-60 60]);                   % -------------------
BL_ShadowPlotBW(ZZb,ZZ(:,[1 2]),DD,lg1,SizeXY,LS([2 1]))
if STAMPA==1; print('-depsc','-painters','-r600',[nomefile 'OG']); end      % -------------------
title('Output Gap - Level','fontweight','bold','fontsize',16)               % -------------------

%%% ===================== %%%
%%%  Figure 5 right plot  %%%
%%% ===================== %%%
ZZ=[cboc(start2:end) chic(:,1)]; ZZ=ML_diff(ZZ(j0:end,:),1);                % Output Gap - 4Q
ZZb{1}=BL_Band(ZZ(:,2),squeeze(ML_diff(chicB(j0:end,1,:),1)),alpha1);       % ---------------
ZZb{2}=BL_Band(ZZ(:,2),squeeze(ML_diff(chicB(j0:end,1,:),1)),alpha2);       % ---------------
SizeXY=ML_SizeXY(Dates4q(j0:end),ML_minmax([ZZ ZZb{1}]),1,5,[-60 60]);      % ---------------
BL_ShadowPlotBW(ZZb,ZZ(:,[1 2]),Dates4q(j0:end),lg1,SizeXY,LS([2 1]))
if STAMPA==1; print('-depsc','-painters','-r600',[nomefile 'OG4Q']); end    % ---------------
title('Output Gap - 4Q % changes','fontweight','bold','fontsize',16)        % ---------------

%%% ========================== %%%
%%%  Figure 4 upper-left plot  %%%
%%% ========================== %%%
YY=[Y2(:,1) cbot(start2:end) chit(:,1)]; YY=YY(j0:end,:);                       % Potential output (levels)
YYb{1}=BL_Band(YY(:,3),(RESTAT.potB(j0:end,:)),alpha1);                         % ------------------------
YYb{2}=BL_Band(YY(:,3),(RESTAT.potB(j0:end,:)),alpha2);                         % ------------------------
SizeXY=ML_SizeXY(DD(1:60),ML_minmax([YY(1:60,:) YYb{1}(1:60,:)]),1,5);          % 70-84
YYY{1}=YYb{1}(1:60,:); YYY{2}=YYb{2}(1:60,:);                                   % ------------------------
BL_ShadowPlotBW(YYY,YY(1:60,:),DD(1:60),lg2,SizeXY,LS([ 3 2 1]),2,[],10)
if STAMPA==1; print('-depsc','-painters','-r600',[nomefile 'PO_7084']); end  	% ------------------------

%%% =========================== %%%
%%%  Figure 4 upper-right plot  %%%
%%% =========================== %%%
SizeXY=ML_SizeXY(DD(61:120),ML_minmax([YY(61:120,:) YYb{1}(61:120,:)]),1,5);    % 85-00
YYY{1}=YYb{1}(61:120,:); YYY{2}=YYb{2}(61:120,:);                               % ------------------------
BL_ShadowPlotBW(YYY,YY(61:120,:),DD(61:120),lg2,SizeXY,LS([ 3 2 1]),2,[],10)
if STAMPA==1; print('-depsc','-painters','-r600',[nomefile 'PO_8500']); end     % ------------------------

%%% ========================== %%%
%%%  Figure 4 lower-left plot  %%%
%%% ========================== %%%
SizeXY=ML_SizeXY(DD(121:end),ML_minmax([YY(121:end,:) YYb{1}(121:end,:)]),1,5);	% 01-18
YYY{1}=YYb{1}(121:end,:); YYY{2}=YYb{2}(121:end,:);                             % ------------------------
BL_ShadowPlotBW(YYY,YY(121:end,:),DD(121:end),lg2,SizeXY,LS([ 3 2 1]),2,[],10)
if STAMPA==1; print('-depsc','-painters','-r600',[nomefile 'PO_0118']); end     % ------------------------

%%% =========================== %%%
%%%  Figure 4 lower-right plot  %%%
%%% =========================== %%%
ZZ=[Y2(:,1) cbot(start2:end) chit(:,1)]; ZZ=ML_diff(ZZ(j0:end,:),1);        % Potential Output - 4Q
ZZb{1}=BL_Band(ZZ(:,3),squeeze(ML_diff(chitB(j0:end,1,:),1)),alpha1);       % ---------------------
ZZb{2}=BL_Band(ZZ(:,3),squeeze(ML_diff(chitB(j0:end,1,:),1)),alpha2);       % ---------------------
SizeXY=ML_SizeXY(Dates4q(j0:end),ML_minmax([ZZ ZZb{1}]),1,5,[-60 60]);      % ---------------------
BL_ShadowPlotBW(ZZb,ZZ,Dates4q(j0:end),lg2,SizeXY,LS([ 3 2 1]))
if STAMPA==1; print('-depsc','-painters','-r600',[nomefile 'PO']); end      % ---------------------
title('Potential Output - 4Q % changes','fontweight','bold','fontsize',16) 	% ---------------------


%%% ===================== %%%
%%%  Figure 1 right plot  %%%
%%% ===================== %%%
id=75;
ZZ=[cboug(start2:end) chic(:,id)]; ZZ=ZZ(j0:end,:); 						% Unemployment gap
ZZb{1}=BL_Band(ZZ(:,2),squeeze(chicB(j0:end,id,:)),alpha1);                 % ----------------
ZZb{2}=BL_Band(ZZ(:,2),squeeze(chicB(j0:end,id,:)),alpha2);                 % ----------------
SizeXY=ML_SizeXY(DD,ML_minmax([ZZ ZZb{1}]),1,5,[-60 60]);                   % ----------------
BL_ShadowPlotBW(ZZb,ZZ,DD,lg1,SizeXY,LS([2 1]))
if STAMPA==1; print('-depsc','-painters','-r600',[nomefile 'UG']); end      % ----------------
title('Unemployment Rate Gap','fontweight','bold','fontsize',16)            % ----------------


%%% ========== %%%
%%%  Figure 7  %%%
%%% ========== %%%
j90=find(year(Dates2)==1990,1);           
ZZ=ML_diff([chit(:,1) chic(:,1) zeta(:,1)],1); YY=ML_diff(Y2(:,1),1);               % Contribution to GDP growth (BL)
ML_ContributionGraphBW(YY(j90-4:4:end),ZZ(j90-4:4:end,:),Dates4q(j90-4:4:end),...     %
    {'GDP','Potential output','Output gap','Idiosyncratic'});                       % --------------------------
if STAMPA==1; print('-depsc','-painters','-r600',[nomefile 'DecompGDP_BL']); end    % --------------------------

ZZ2=ML_diff([cbot(start2:end,1) cboc(start2:end,1)],1); YY=ML_diff(Y2(:,1),1);      % Contribution to GDP growth (CBO)
ML_ContributionGraphBW(YY(j90-4:4:end),ZZ2(j90-4:4:end,:),Dates4q(j90-4:4:end),...    %
    {'GDP','Potential output','Output gap'});                                       % --------------------------
if STAMPA==1; print('-depsc','-painters','-r600',[nomefile 'DecompGDP_CBO']); end   % --------------------------

